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ABSTRACT 

We present results from large-scale particle simulations of the viscous overstability in Saturn's 
rings. The overstability generates a variety of structure on scales covering a few hundred 
metres to several kilometres, including axisymmetric wavetrains and larger-scale modulations. 
Such patterns have been observed in Saturn's rings by the Cassini spacecraft. Our simulations 
model the collisional evolution of particles in a co-rotating patch of the disk. These are the 
largest A^-body simulations of the viscous overstability yet performed. The radial box size is 
five orders of magnitude larger than a typical particle radius, and so describes a 20-50 km 
radial portion of the rings. Its evolution is tracked for more than 10000 orbits. In agreement 
with hydrodynamics, our A^-body simulations reveal that the viscous overstability exhibits a 
rich set of dynamics characterised by nonlinear travelling waves with wavelengths of a few 
hundred meters. In addition, wave defects, such as sources and shocks, punctuate this bed 
of waves and break them up into large-scale divisions of radial width ~ 5 km. We find that 
the wavelength of the travelling waves is positively correlated with the mean optical depth. 
In order to assess the role of the numerical boundary conditions and also background ring 
structure, we include simulations of broad spreading rings and simulations with a gradient in 
the background surface density. Overall, our numerical results and approach provide a tool 
with which to interpret Cassini occultation observations of microstructure in Saturn's rings. 
We present an example of such a synthetic occultation observation and discuss what features 
to expect. We also make the entire source code freely available. 
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> 1 INTRODUCTION 
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Saturn's rings exhibit a phenomenal wealth of axisymmetric struc- 
ture that covers a vast range of lengthscales, from 100 m to 100 km. 
On the shortest scales, ultraviolet and radio occultation observa- 
tions reveal that the variations are quasi-periodic with wavelengths 
between 150 m and 220 m. These wavetrains intermittently popu- 
late regions of optical depth ~ 1-2, such as in the A-ring and less 
opaque areas in the B-ring ( |Colwell et al.|[2007] |Thomson et al.| 
|2007[ [Colwell et al.|2009| l. In addition, the Cassini cameras have 
uncovered disjunct bands of varying brightness with widths of 1-10 
km. These intermediate-scale variations often coincide with the pe- 
riodic microstructure (Porco et al. 2005 ). On the other hand, broad 
undulations on 100 km scales pervade the inner B-ring, while the 
central B-ring splits into alternating zones of low and high optical 
depth, also on scales of 100 km (Col well et al.|20 09 ). We seek to 
understand how and why these different features develop and the 
nature of any relationship between them, while concentrating on 
the structure on short and intermediate scales. 



Fine-scale structure originates in a linear instability connected 
to the viscous stress of the ring, termed the 'viscous overstability'. 
Its fastest growing modes take the form of axisymmetric density 
waves on wavelengths ~ 50 m (Schmit & Tscharnuter 1995 1. Such 
waves give rise to perturbations in the viscous stress that extract en- 
ergy from the background orbital shear. This surplus energy is fed 
into the wave, where it outcompetes standard viscous damping, and 
thus leads to runaway growth. We concentrate in this paper on the 
small-scale axisymmetric manifestation of the instability, but note 
that the same mechanism can lead to growth of a global eccentric 
mode in both narrow and broad rings (Borderies et al. 1985] |Pa 



I). 



|paloizou&Lin[T 988 Long aretti & Rappaport|1995[|Ogilvie|2001 

The viscous overstability has been studied in the context of 
continuum models (both hydrodynamical and kinetic) and Af-body 
simulations, with the early stages of its evolution receiving the 
greatest attention ( Schmit & Tscharnut er|1995[|Salo|2001||Schmidt| 
|et al.|2001||Salo et al.|2001||Latter & Ogilvie|2006a| |2008) . The 
hydrodynamical studies reveal that instability occurs once the ring 
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viscosity increases sufficiently steeply with surface density; i.e. if 
P = dlnv/dlnr > 1, where v is viscosity and r is normal opti- 
cal depth. In kinetic models and A'-body simulations this instability 
condition corresponds to r > r c , where t c is a critical optical depth. 
The latter theoretical formalisms, however, yield no upper bound on 
t, above which overstability disappears, as is indicated by Cassini 
observations of the B-ring. The nonlinear saturation of the insta- 
bility, on the other hand, has mainly been explored with hydro- 
dyna mical models ([ Schmit & Tscharnuter 1999; Sch midt & Salo| 
|2003||Latter & Ogilvie|2009||2010| >. The most recent studies show 
that the long-term evolution of the ring is dominated by nonlinear 
travelling wavetrains with wavelengths ~ 200 m — shorter waves 
are vulnerable to secondary instabilities. The simulated wavetrains 
undergo small chaotic fluctuations and can be punctuated by wave 
defects (sources and shocks) that partition the radial domain into 
structures on larger scales > 1 km (Latter & Og ilvie 2010}. The 
similarity between the shorter wave features and the larger scale 
divisions, on the one hand, and the diverse Cassini observations, 
on the other, is encouraging, though additional work is necessary 
to make the correspondence tighter. In particular, the role of self- 
gravity has yet to be clarified. 

In this paper we simulate the nonlinear saturation of the vis- 
cous overstability with large-scale particle simulations in a local 
co-rotating patch of disk (the shearing box approximation). Our 
simulations generalise the previous hydrodynamical results to the 
more realistic regime of a dense granular flow, in which the ring's 
pressure tensor deviates significantly from the Newtonian prescrip- 
tion (Latter &~O gilvie 2006a, 2008 ). They simultaneously extend 
previous A'-body simulations to the greater length and time-scale 
necessary to capture the nonlinear dynamics properly. These scales 
are now computationally feasible thanks to the freely available col- 
lisional A'-body REBOUND (Rein & Liu 2012) which includes a sym- 
plectic integrator for the shearing sheet and an efficient collision de- 
tection algorithm. We omit self-gravity; its influence will be tested 
in future work. A beneficial consequence of this omission is that 
we can shorten the azimuthal extent of the box. This permits us to 
simulate boxes with radial sizes of 55 km (55 000 particle radii) 
for some 10000 orbits, orders of magnitudes more than previous 
simulations. 

Our results are in general qualitative agreement with the hy- 
drodynamical simulations of |Latter & OgilvTe] (2010) . Both ap- 
proaches yield simulations that are dominated by travelling non- 
linear wavetrains of comparable morphology and wavelength (~ 
100 m), though the particle simulations exhibit a protracted early 
stage in the evolution in which the counterpropagating waves inter- 
act in a complicated way, giving rise to standing wave figures and 
beating patterns ('wave turbulence'). In addition, the A'-body sim- 
ulations yield sources and shocks structures as in hydrodynamics 
though their sizes and morphology differ. To investigate the poten- 
tial role of global ring structure (and to minimise that of the pe- 
riodic boundary conditions), we also performed simulations with 
viscously spreading rings (mimicking open boundaries) and rings 
with a background density gradient. We stress that the structures 
that emerge in our simulations never exhibit lengthscales longer 
than roughly 5-10 km. We conclude that the observed 100 km fea- 
tures in the B-ring are unlikely to be generated by viscous oversta- 
bility. 

New Cassini observations of star occultations have a resolu- 
tion that is below a few hundred meters such that individual peaks 
of the overstability-supported wavetrains can be directly observed 



(M. Hedman, private communication). A direct comparison to nu- 
merical simulations could allow for the measurements of particles 
properties such as size, density and coefficient of restitution. We 
present the first step towards this goal and present synthetic oc- 
cultation observations from our A'-body simulations that show the 
feasibility of this approach. 

The paper is organised as follows. In Section[2]we describe the 
physical model and governing equations, and then the numerical 
method we use to solve them. Section|3]provides a brief comparison 
between our code and the kinetic equilibria and linear theory of 
|Latter & Ogiiyie] p008] >, as well as a number of numerical tests. 
In Section[4]we present our main results, beginning with a detailed 
analysis of a large-scale fiducial simulation and then moving on 
to the role of the main parameters and global disk structure. We 
summarise our results and discuss the implications in Section]?] 



2 PHYSICAL AND NUMERICAL MODEL 
2.1 Equations of motion 

We solve the equations of motion in the Hill approximation (Hill 
|1878[ > which is a local coordinate system that is co-rotating with a 
particle on a circular orbit. The gravity from the central object is 
linearised in local co-ordinates and the orbital frequency is a con- 
stant. This allows, but does not restrict, us to shear-periodic bound- 
ary conditions. In that case the Hill approximation is also referred 
to as the shearing sheet. The x, y, and z coordinates point in the 
radial, azimuthal, and vertical direction respectively. The equations 
of motion for a test particle can then be written as 

x = lily + 3Q. 2 x 

y = -2Slx (1) 
z = -Orz, 

where Q and SI- are the angular velocity and vertical epicyclic fre- 
quency, respectively. The solution to these equations can be written 
as epicycles (e.g. |Rein & T remaine 2011). In the case where the 
central object is a point mass and self-gravity is neglected, we have 
£2 = fl z . We denote the radial length of the box by L x and the 
azimuthal length by L y . Typically L y <k L x in the simulations we 
perform. 

The only further ingredients needed besides Eqs. |TJ are the 
finite particle size r p and a collision model. We treat particles as 
hard spheres (they are not permitted to overlap) and the outcome 
of a collision is described using a normal coefficient of restitution 
6, which can be either constant or a function of the impact velocity 
Vj mp . The particles have no spin. 

In this paper we do not treat self-gravity directly. This is done 
only for simplicity and clarity. However, we increase the vertical 
epicyclic frequency in some simulations Sl z = 3.6Q in order to con- 
centrate particles in the mid-plane, thus mimicking self-gravity's 
vertical compression. This convenient method has been used in 
most previous ring studies (e.g. .Wisdom & Tremaine 1988] [Salo] 
|et al.|200T) . It allows us to remove one scale from the problem (i.e. 
the particle density or Hill radius) and more easily focus on the 
fundamental details of the collective axisymmetric dynamics. 

We plan to study simulations including self-gravity in a fu- 
ture paper. Self-gravity is important because it may instigate non- 
axisymmetric wake structures that could compete with the oversta- 
bility and alter its saturation (Salo et al. 2001). Self-gravity will fur- 
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thermore change the properties of the overstability waves directly 
through their nonlinear dispersion relation. 

2.2 Diagnostics 

In order to probe the collective behaviour of the granular flow, we 
require a number of averaged quantities. We define the mean nor- 
mal geometrical optical depth f as the total projected area of the 
particles on the (x, y) plane divided by the total area of the (x, y) 
plane. In other words 

f = Nnr 2 p l(L x L y ), (2) 

where N is the number of particles. Thus, f is stipulated at the be- 
ginning of each run and will not change. 

We also define the radially and temporally varying optical 
depth, denoted simply by t, by subdividing the radial domain into 
thin strips of radial length L s : 

T(x i ) = N,Tcr 2 p /(L s L y ), (3) 

where x, is the radial location and A', is the number of particles 
in the i'th strip. As we typically begin a run from a homogeneous 
state, at t = we have t(x,) = f for all i. The ensuing overstable 
fluctuations will subsequently be captured by the variations in t(x;). 
L s is chosen to be smaller than any scale we are interested in but 
not so small as to be influenced by Poisson noise. 

We further introduce the photometric optical depth f which we 
define as one minus the transmission. The transmission (the frac- 
tion of light passing through the ring) is defined to be 1 for f = 
and for f = 1. This is calculated with a Monte Carlo Ray tracing 
algorithm. The algorithm shoots rays through a random point in the 
simulation box but with a predefined direction (the viewing angle). 
Then, we check for intersections with ring particles along the ray. 
This is repeated until a desired accuracy has been reached. Unless 
otherwise noted we use a viewing angle normal to the ring plane. 

The filling factor is defined as the proportion of volume taken 
up by the particles. For spherical particles it can be defined as 
FF = (4n/3i)nri, where n is volumetric number density. Particu- 
larly useful is the filling factor at the midplane FF , which requires 
the calculation of the number density at z = 0. 

Let us further introduce the mean velocity dispersion tensor 
which is computed from 

Wij = (x iXj h (4) 

where (xi,X2,X3) = (x, y, z) and the angle brackets indicate a suit- 
able average over the particles and possibly over time. The velocity 
dispersion c 2 is then Wa/3. The translational (local) component of 
the kinematic viscosity is simply 

The collisional (nonlocal) component of the viscosity is 

Vco " = 3flMT 2 (x > _ *< )Ap " (6) 
where the sum is taken over all binary collisions that occur in a 
time interval T. Here M is the total mass of all ring particles, Ap y 
denotes the transfer of y momentum from the inner to the outer 
particle in such a collision, and x> and x< the radial locations of the 
two impacting particles (Wisdom & Tremaine 1988 , Daisaka et al. 
|2001[ >. As we neglect self-gravity there is no gravitational or wake 
contribution to the overall momentum transport. 



2.3 N-body code and algorithms used 

We use the freely available A'-body code REBOUND ( Rein & Liu| 
|2012| > to perform all of the simulations presented in the paper. The 
code is ideally suited for this kind of study because of both speed 
and accuracy. We use the highly efficient and accurate mixed vari- 
able symplectic integrator SEI (Rein & Tremaine 201 1 ). Collisions 
are detected using a so called plane-sweep algorithm. The plane- 
sweep algorithm scales linearly with the number of particles, O(N), 
for highly elongated boxes such as those considered here. Both al- 
gorithms are already implemented in REBOUND. These advantages 
allow us to use boxes which are 10 5 times larger than the typical 
particle radius. We therefore do not need to rely on periodic bound- 
ary conditions and can study effects near ring boundaries and den- 
sity gradients. 

One crucial aspect in simulating the viscous overstability cor- 
rectly is to make sure that there is no preferred direction introduced 
by the numerical scheme. In fact, the symplectic integrator intro- 
duced by |Rein & Trem aine (2011) ensures this symmetry is not 
broken, even on the scale where machine precision becomes im- 
portant. However, collisions formally break the symplectic nature 
in our implementation. Therefore, we randomise the order in which 
collisions are resolved after each timestep. This has proven to be 
the cleanest and most efficient way to remove all spurious corre- 
lations which might otherwise be introduced in the order in which 
collisions are resolved. Unfortunately, this complicates the paral- 
lelisation on distributed memory systems. 

For the reader's convenience, we list all the simulations pre- 
sented in this paper in Table [T] The first column describes the sim- 
ulation. The second column lists the particle radius. The third and 
fourth columns list the radial and azimuthal size of the box. The 
fifth and sixth columns list the initial mean optical depth f and the 
number of particles. The seventh column lists the vertical epicyclic 
frequency in units of the standard epicyclic frequency. The eighth 
column lists the coefficient of restitution. 



3 COMPARISONS AND NUMERICAL TESTS 
3.1 Background equilibrium state 

We begin by calculating the properties of homogenous statistical 
equilibria with our numerical code and comparing them to the ki- 
netic theory of Araki & Tremaine ( 1986) as reformulated by [Latter] 
|& Ogilvfe] < |2008) >. We use a simulation box with a width of only 
25 particle radii. As a consequence, growing overstability modes 
cannot fit into the domain (see also Sect. \3.3\ . Depending on the 
optical depth, the number of particles ranges from N = 20 to 600. 
A constant coefficient of restitution e = 0.5 is used. Other than the 
particle radius r p and the frequencies CI and Q.,, which we set equal 
to 1 in our units, there is no other scale. 

In Figure [T] we plot the filling factor in the midplane FF , 
the velocity dispersion c 2 , the xx, yy and xy components of the 
velocity dispersion tensor W,-j, and the total kinematic viscosity 
v totai = v coii + y trans- These quantities are plotted as functions of the 
optical depth f . Overall, the numerical equilibria are in good agree- 
ment with the kinetic theory, with errors of a similar order to pre- 
vious comparisons (Wisdom & Tremaine 1988, Latter & Ogilvie 
|2008| >. At larger f , systematic discrepancies emerge which are the 
result of the large filling factors FFq x 0.4. In this regime the En- 
skog approximation begins to break down, and the kinetic results 
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Name / Description 
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3.6 
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15.53 


1.64 
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3.6 


0.5 
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18641 


15.53 


1.64 
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3.6 
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Bridges et al. 
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1.64 
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Table 1. Parameters of simulations presented in this paper. See text for details. 




0.5 1 1.5 2 2.5 0.5 1 1.5 2 2.5 3 0.5 1 1 .5 2 2.5 3 

optical depth t optical depth % optical depth t 




optical depth t optical depth z optical depth t 

Figure 1. Comparison of numerical simulations to kinetic theory. From top left to bottom right: filling factor, total viscosity, velocity dispersion, xx, yy and xy 
components of the velocity dispersion tensor. See text for details. 



are less reliable. There is also a 10% error in the xy component 
of the viscous stress at lower f , which has been previously noted 
( |Latter &T3 gilvie 2008 ). These results show that our particle code 
correctly reproduces established kinetic ring equilibria. 



3.2 Linear stability 

Next we examine how well the code treats the linear stability of 
the equilibrium states. Instead of examining the growth rates of the 
viscous overstability modes, as in Schmidt et al. (2001), we de- 
scribe the curve of marginal stability in the full parameter space of 
6 and f . This result extends previous W-body work that tested only 
a handful of cases l |Salo et al.|200T] >, and can directly connect to the 
predictions of the kinetic modelling (Latter & Ogilvie 20081. 

For given collisional properties, it has been shown that vis- 



cous overstability emerges once the optical depth increases above 
a critical value. This f dependence is actually a proxy for central 
filling factor FFq, which more directly influences the behaviour of 
the ring's viscous stress. As FF increases, so does the sensitivity 
of the stress to density variations, as mediated by the Enskog fac- 
tor in the kinetic theory ( Araki & Tremaine 1986; Latter & Ogilvie 
2008). If the stress increases sufficiently fast with density (i.e. if j5 
is sufficiently large) then instability ensues ( Schmit & Tsch arnuter] 
[7995]|Schmidt et al.|2001> . 

We compute the numerical stability criterion in the following 
way. A grid of f and e values is set up and a sequence of simu- 
lations run, each corresponding to a unique (f, e). The simulations 
use L x = 1000 m and are evolved for up to 10 4 orbits, which en- 
sures that there is sufficient time for the unstable modes in the box 
to achieve appreciable amplitudes. The level of activity in the box 
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0.5 1 1.5 2 2.5 

optical depth t 



Figure 2. Deviation from the initial optical depth as a function of the initial 
geometric optical depth and the coefficient of restitution. This is a proxy for 
the activity of the viscous overstability. 

is measured by ((r(f) - f) 2 ), where the angle brackets denote a box 
average. In order to obtain sufficiently large FF for instability at 
a reasonable f we amplify fl ; to 3.6 CI. Doing so mimics the com- 
pression of the disk via self-gravitation (see above). The result of 
these simulations are plotted in Fig. [2] in which we show the mag- 
nitude of the activity at the end of the run. Superimposed on the 
figure in red is the stability curve calculated from the dense gas 
kinetic theory ( |Latter & Ogilvie|2008"l >. 

As any activity corresponds to the onset of viscous overstabil- 
ity, the purely black regions in Fig.[2]signify stability. All non-black 
regions indicate instability, except for a handful of outliers influ- 
enced by noise. Clearly, for given e less than roughly 0.9 there is a 
critical f above which we have instability, in agreement with previ- 
ous results. Larger e generally requires a larger f . This is because 
the velocity dispersion and disk thickness tend to increase with e, 
there being less dissipation. As a consequence, greater filling fac- 
tors are more difficult to achieve at lower optical depths. There also 
appears to be a hard cut-off around f = 0.5 below which instability 
can never occur. 

There is general qualitative agreement with the kinetic 
marginal curve. The quantitative discrepancies arise because insta- 
bility occurs for larger FFo (above roughly 0.35), the regime in 
which the kinetic theory is less reliable. For larger e this is espe- 
cially the case, as here FF approaches 0.5 so as to force instability 
in these less dissipative conditions. 

In summary, both the equilibrium and linear stability results 
show that our particle code is behaving in the correct manner, with 
regards to previous calculations. This gives us confidence in its ca- 
pability to describing the nonlinear evolution of the viscous over- 
stability which is the main focus of this paper. 

3.3 Numerical convergence 

We verify how the results of the simulations depend on the size 
of the box and the timestep. These are the only two numerical pa- 
rameters we need to check. In these simulations the geometric op- 
tical depth is f = 1.64 and we use a constant coefficient of resti- 
tution e = 0.5. As we have shown in Sect. |3.2| this system will be 
viscously overstable. 

We integrate the system for approximately 320 orbits which 



gives the overstability enough time to develop. The power spec- 
tra of the radial variations in the density, the average height of a 
particle above the mid-plane and the average velocity of a particle 
are then measured 100 times per orbit and averaged. One of these 
quantities, the power spectrum of the density, is plotted in Fig. [3] 
for different dt, L x and L y . The other spectra are qualitatively very 
similar. 

In Fig. |3(a)| we fix L x = 1024m and L y = 5m and vary the 
timestep. One can see that a timestep dt = 2n ■ lO ^ST 1 or smaller 
is sufficient to capture the overstability accurately. From Fig. |3(b)| 
where all parameters except the azimuthal box width are fixed, it 
is clear that L y < 4m suppresses the overstability. For L y ^ 10m 
the spectra are independent of Ly. Thus the box size is sufficiently 
large to capture the overstability accurately. Finally, by inspecting 
Fig. |3(c)| one can see that a radial box width of L x = 1024 m is suf- 
ficient to accurately capture the dominant wavelengths in the non- 
linear saturated state. 

To summarise, the development of the viscous overstability in 
these smaller runs is relatively insensitive to the numerical param- 
eters. Of course this does not tell us whether the large-scale nonlin- 
ear behaviour is captured adequately, but it nonetheless provides a 
foundation of confidence that this is indeed the case. 

We should also point out that in simulations with higher opti- 
cal depth, f ^ 2, we need a larger L y to avoid crystallisation of par- 
ticles. Quasi-crystalline structures are a lower energy state and act 
as an attractor if the velocity dispersion is too small to randomise 
the particle motion. Higher density rings are more prone to this ef- 
fect because a higher collision rate leads to a cooler and denser ring 
with strong excluded volume effects and particle correlations. This 
is purely an artefact of using equal-sized particles and will not oc- 
cur in a real ring£] We could simply allow for a size distribution to 
avoid this problem (we present one such simulation in Sect. |4.2| . 
However, this would complicate the analysis and comparison to an- 
alytical estimates. 



4 NUMERICAL RESULTS 

In this section we present our main results on the nonlinear satu- 
ration of the viscous overstability in large radial domains and over 
long times. We ran a significant number of simulations varying our 
two main control parameters f and e, but we will describe only one 
set of parameters in detail. This 'fiducial simulation' conveniently 
exhibits all the main features witnessed in the other runs. In addi- 
tion, we tested how different global structures impacted on the long 
term behaviour. Specifically, simulations are shown of a spreading 
ring and a ring with a density gradient. 



4.1 Fiducial simulation 

Our illustrative run employs an optical depth of f = 1.64. We use 
shearing sheet boundary conditions with an azimuthal and radial 



Analogous crystallisation phenomena have been hypothesised as a re- 
sult of the (poorly constrained) surface adhesion between ring particles 
by |Trem aine (20031. See also |Perrine et al.|(201 l) ; |Perrine & Richardson| 
|2012} . 
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(a) Density power spectra for different timesteps 
dt. The unit of the timestep is 2;rfl~' . 



(b) Density power spectra for different azimuthal 
box widths L v . 



(c) Density power spectra for different radial box 
widths L, . 



Figure 3. Convergence study. Power spectra of the density for different timesteps and box sizes. 




1000 2000 3000 4000 5000 6000 7000 8000 
time [orb] 

Figure 4. Kinetic energy evolution (only x and y components) in the fiducial 
simulation. 



width of Ly = 15.53 m and L x = 12.43 km, respectively. All parti- 
cles are initially randomly distributed. The simulation was then in- 
tegrated for 8200 orbits which corresponds to over 10 Earth-years, 
assuming a distance from Saturn of a = 130000 km. 

Viscous overstability attains nonlinear amplitudes after a few 
hundred orbits, and then settles into a disordered state for some 
5000 orbits. After roughly 5300 orbits this state resolves into a sin- 
gle travelling wavetrain of relatively uniform amplitude. The first 
disordered stage of the evolution can be arranged into two con- 
secutive phases, (a) a 'wave turbulence', characterised by nonlin- 
early interacting travelling waves, standing waves, and beating pat- 
terns, and (b) a cellular source/shock state, in which relatively well- 
defined patches of counterpropagating waves are separated by wave 
defects (shocks and sources). In Fig.|4]we present the total fluctu- 
ating kinetic energy of the box as a function of time. In it the three 
phases of the evolution can be distinguished. From roughly 500 
orbits to 2000 orbits, there is a linear increase of energy associ- 
ated with stage (a), the 'wave turbulence'. Then from orbits 2000 
to 5300 the energy saturates, while undergoing long-term oscilla- 
tions. This corresponds to phase (b), the source/shock state. Af- 
ter 5300 orbits the system relaxes to a single travelling wavetrain 
that fills up the entire box and yields a single uniform energy. This 
is the end-state of the evolution and is shared by all simulations 
that used shear-periodic radial boundaries. We now discuss each of 
these three phases in more detail. 

We begin in reverse with the final state, which we plot in 
Fig. |5(a)| This figure shows a snapshot of the particle positions in 
the (x,y) plane and a (x, z) cross-section, the optical depth t, the y- 
and z-averaged azimuthal velocities (relative to the shear) v y , and 



the velocity dispersion c after / = 8000 orbits. We find that the fi- 
nal stable nonlinear wavetrain state propagates in a single direction. 
The first plot confirms that the nonlinear evolution remains axisym- 
metric (within the limited azimuthal extent of the box), while the 
r profile exhibits the characteristic cusp-like shape shared by all 
density waves in thin disks, where particle motions are controlled 
by epicycles. The system has yet to completely relax to a wave- 
train of uniform amplitude, and displays a large-scale modulation 
in wavelength. Ultimately this modulation decays. Note that the 
wavelength variation corresponds to an amplitude variation in the 
velocity field, as shown by the hydrodynamical theory. The aver- 
age wavelength exhibited is some 450 m. The velocity dispersion 
profile reveals spikes in the random motions at wavetrain peaks. 
It is in these high density peaks that the orbital energy is primar- 
ily thermalised, as the collision frequency (and associated colli- 
sional stresses) are maximised there. The (x, z) cross-section shows 
these regions of high-density are accompanied by coherent vertical 
expansions (splashing) via excluded volume effects. This vertical 
structure is connected to the vertical 'breathing mode' exhibited by 
viscous overstability in gaseous disks {Latter & O gilvie 2006b]!. 

The source/shock state is described in Fig. |5(b)| which is sim- 
ilar to the previous figure but omits the v v and c profiles. The t 
profile shows the existence of two patches of density waves of vary- 
ing wavelength with their direction of propagation indicated by ar- 
rows. The two wavetrains are separated by a source region located 
at roughly x = m, which generates waves, and a shock region lo- 
cated around x = -9000 m, which is where the counterpropagating 
waves collide and destroy each other. The shock region straddles 
the (periodic) radial boundary. This configuration of dynamical el- 
ements reflects behaviour in the complex Ginzburg-Landau equa- 
tion, which serves as a good reduced model for the overstable sys- 
tem {Aranson & Kramer|2002"l |Latter & Ogilvie|200"9"l >. It also re- 
produces qualitatively the hydrodynamical simulations of Latter & 
Ogilvie (2010), yet there exist interesting discrepancies. The wave- 
trains themselves are far more disordered and undergo larger am- 
plitude fluctuations in their amplitude and phase. In addition, the 
source and shock regions are broader and less defined. In particular, 
the shock region can extend over 1000-1500 km, i.e. many char- 
acteristic wavelengths of the wavetrain. This is because colliding 
wavetrains penetrate each other more successfully in the particle 
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(a) Fiducial simulation after t = 8000 orbits. 
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(b) Fiducial simulation after t = 4000 orbits. 




(c) Fiducial simulation after t = 300 orbits. 
Figure 5. Snapshot of the fiducial simulation at different times. 
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simulations and must undergo multiple 'collisions' before they dis- 
sipate. By contrast, the hydro simulations support shocks that are 
very narrow and well-defined: waves are almost totally destroyed 
in a single collision. This difference in the nonlinear interaction 
between travelling wavetrains is perhaps the key disagreement be- 
tween the hydrodynamical and A'-body modelling. 

In Figure [5(c)! we plot a snapshot of the evolution during the 
initial disordered wave state. The dynamics in this phase is compli- 
cated, comprising stretches of counterpropagating travelling waves 
interpenetrating each other and giving rise to various nonlinear in- 
teractions such as beating and standing wave patterns. It is possi- 
ble that source and shock regions exist in the wave maelstrom but 
they are so radially extended that they overlap significantly and do 
not function in the same way as in the shock/source stage. As is 
clear by comparison with Fig. |5(b)| characteristic wavelengths are 
much shorter in this evolutionary stage, closer to the wavelength 
of the most unstable mode (~ 100 m), and they steadily increase 
with time ( Schmit & Tscharnuter 1999 ). This, in fact, drives the in- 
crease in kinetic energy witnessed in Fig. [4] as the wave amplitude 
is correlated with wavelength.. This stage of the evolution bears 
similarities with the one-dimensional 'wave turbulence' that char- 
acterises certain plasmas, nonlinear optics, and other physical phe- 
nomena (Nazarenko 20TTj for e.g.). Notably, the overstable ring 
shares with those systems an inverse cascade of energy to longer 
lengths and a general movement towards greater nonlinearity and 
coherent structure. 

The stages of the dynamical evolution can be summarised in 
a stroboscopic space-time diagram, Fig. [6] in which the optical 
depth is plotted as a function of x and f, but only at every 10th 
orbit. This screens out the fast time associated with the wave phase 
speed yet keeps the sense of their and the group velocity's direc- 
tion, while also indicating important features such as the sources 
and shocks. The first clear signal of a source region emerges at 
around t = 1800 orbits near x = -5000 m. It can be distinguished 
by the light diagonal rays propagating away on either side. These 
rays represent the slower group velocity of the waves (not the very 
rapid phase velocity ~ 100m per orbit). The shock region is less 
easy to pick out. At r = 2000 orbits it extends from x = - 100 m to 
x = 1000 m. As illustrated by hydrodynamical simulations, source 
and shock regions migrate around the domain on a slow time-scale. 
Eventually they collide and annihilate each other. This occurs near 
t = 5300 orbits. After this point we have a single wavetrain prop- 
agating in one direction (all the rays point in the same direction). 
Another interesting feature is that the wavelength of the waves gen- 
erated by the source vary with time and can be different depending 
on which side of the source they originate. The time-variable be- 
haviour was noted in the hydrodynamical simulations but was less 
prominent; the source asymmetry appears a novel feature of the 
A'-body simulations. 

In all the runs undertaken with shearing periodic boundaries 
the long term saturation of overstability is a single travelling wave- 
train of uniform amplitude. However, this final state is to some 
some extent an artefact of the shearing box's translational sym- 
metry in radius. As far as the collective axisymmetric dynamics 
are concerned, there is no preferred radius. A real ring, on the 
other hand, exhibits radial structure in particle size, composition, 
and number density on intermediate to long scales. Moreover, res- 
onances and wakes create spatially fixed perturbations. These vari- 
ations should prohibit the transition to a smooth and coherent state 
with one wave-train. It is therefore likely that the real rings will ex- 




Figure 6. Stroboscopic space-time plot of the fiducial simulation. 



hibit the dynamics exhibited by simulations on intermediate times, 
before the influence of the boundary conditions dominate — in 
other words, before, the box realises that it is periodic. Viscous 
overstability in the real rings will probably manifest as either (a) the 
disordered wave state or (b) the source/shock state. Perhaps differ- 
ent regions in the rings support one or the other state, giving rise to 
differing observations on intermediate scales 1-10 km. As revealed 
by the Cassini cameras, there exists a rich and irregular set of radial 
variations on these scales (Porco et al. 2005 ). In Section|43]we trial 
other global set-ups and boundary conditions in order to break the 
radial translational symmetry and to assess how the overstability 
saturates in the presence of simple global structures. 

We ran one additional simulation with the same parameters as 



© 0000 RAS, MNRAS 000, 000-000 



Large-scale N-body simulations of the viscous overstability in Saturn's rings 9 




100 200 300 400 500 600 
wave length [m] 

Figure 7. Power (arbitrary units) as a function of wavelength (horizontal 
axis) and the optical depth f (vertical axis). Each row corresponds to one 
simulation. 

the fiducial simulation but with an extremely large radial domain, 
spanning over 55 000 particle radii (55 km). More details of the 
larger simulation, as well as a stroboscopic space-time diagram, are 
presented in Appendix [A] As is discussed there, precisely the same 
phenomena appears on similar scales in the large box. This shows 
the box independence of our results. It also suggests that the longest 
scales generated by viscous overstability do not exceed 5-10 km. 

4.2 Exploration of different parameters 

In order to test the robustness of the results presented in the pre- 
vious subsection, we conducted simulations with different e and 
f . In particular, we used both e = (no bouncing) and the ve- 
locity dependant Brid ges et al.| (T984) collision law while keeping 
t = 1.64. Neither resulted in nonlinear evolutions significantly dif- 
ferent to that of the fiducial run. This makes us confident that, for 
sufficiently inelastic collisions, the precise treatment of the coeffi- 
cient of restitution is unimportant and that a value of 6 = 0.5 yields 
representative results. It is likely that larger e may lead to different 
nonlinear behaviour, as the ring will generally be 'hotter' and less 
dense. However, experimental and theoretical evidence points to 
very dissipative collisions ( |Porco et"aL 2008 Sch midt et al.| 2009), 
so this hot regime is probably irrelevant to conditions in Saturn's 
rings. 

We undertook additional large-scale simulations with two dif- 
ferent f, equal to 1 and 2. In each case e = 0.5. In neither case, 
was the general structure of the nonlinear evolution different to the 
fiducial case. It should be noted, however, that the lower f run ex- 
hibited less violent chaotic fluctuations and more ordered wave- 
trains during the source/shock phase. What these simulations also 
revealed was that the characteristic wavelength of the dynamics (in 
which the power was concentrated) positively correlated with opti- 
cal depth. Lower f runs yield shorter wavelength behaviour, while 
higher f gives a longer characteristic wavelength. This correlation 
was also noted indirectly in Latter & Ogilvie ( 2009) in their Table 
3, which shows that the wavelength of the first stable wavetrain in- 
creased with f (as mediated by the dimensionless parameter P). It 
is likely that in the disordered nonlinear state, this first stable wave- 
length works as an attractor and fixes the characteristic wavelength 
on sufficiently long times. 



To develop this idea we subsequently ran a sequence of 
smaller scale simulations for different f, in which L, = 4000 m. 
Each simulation was run for 5000 orbits. Because of the smaller 
box sizes, after this length of time the system has resolved into a 
single stable wavetrain. Fig.|7]summarises these results by plotting 
the power in each wavelength for each f . As is clear, the wave- 
length of most power increases linearly with optical depth, begin- 
ning at 210 m at f = 1 and rising to roughly 400 m at f = 2. Note 
that the viscous overstability is not active for f < 0.75, consistent 
with the results in Sect. |3.2| A correlation of wavelength with opti- 
cal depth might have been observed in recent Cassini observations 
(M. Hedman, private communication). 

All simulations presented so far use a single size of ring par- 
ticles. We also ran an additional simulation using a size distribu- 
tion. The distribution has a slope of dN/dr p cc r~ 3 and r p ranges 
from 0.5 m to 2 m, motivated by current estimates of the typical 
size distribution in Saturn's rings ( Cuzzi et al. 2009 1. The system is 
viscously overstable and the non-linear behaviour is indistinguish- 
able from the fiducial simulation. As an illustration, we present a 
3D rendering of a small part of the (x, z) plane in Fig. [8] The plot 
shows approximately two wavelengths. Clearly visible is the verti- 
cal splashing of particles. These simulations make us confident that 
our results are general and are not an artefact of the single sized 
ring particles. There might be quantitative differences that depend 
on the precise nature of the size distribution. However, our results 
indicate that these effects will be small and we therefore do not 
investigate them any further in this paper. 

4.3 Exploration of different global structures 

In this subsection we extend the shearing box model to account for 
different global structures: a large spreading ring and a ring with 
a radial density gradient. We then check how these influenced the 
long term saturation of the overstability. The large extent of our 
simulations (covering scales separated by almost 5 orders of mag- 
nitude — 1 m to 50 km) allows a convincing exploration of global 
structure in direct A'-body simulations. 

4.3.1 Open boundary conditions 

By open boundary conditions we mean a ring that has a finite ra- 
dial width less than the radial box size L x . Thus the ring can spread 
viscously. The boundaries are, in fact, not really open: if particles 
exit the box on one side, they enter the box on the other side. How- 
ever, we initially place the box boundaries far from the disk edge 
(2000 m) and stop the simulation when the ring has spread to the 
point that particles cross the boundary. 

We present a space-time plot of such a simulation in Fig- 
ure [9(a)] One can see that the ring edges are indeed spreading vis- 
cously as expected (i.e. ~ V0- The development of the overstability 
in the middle of the ring is initially unchanged by these boundary 
conditions: the instability remains local on early and intermediate 
times. This is aided by the fact that overstability-generated wave- 
trains are not reflected by the ring edges. Waves, once they reach the 
edge, are lost to the system and so the spreading edge behaves like 
an open boundary or the buffered boundaries employed in |Latter &| 
Ogilvie (2010). Near t = 500 orbits two sources have manifested; 
one at x = -4000 m and the other at x = 4000 m. After 1200 or- 
bits, after a complicated interaction, the intervening shock region 
annihilates itself and the leftmost source leaving only the rightmost 
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Figure 8. 3D rendering that shows a small part of the (x, z) plane in the simulation using a particle size distribution. 
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Figure 9. Stroboscopic space-time diagram of different global structures. 



source. This remaining source stays at this locations for the rest of 
the simulation, as in the hydrodynamical simulations of buffered 
boundaries. 

This persistence of a single source is the key difference to sim- 
ulations using the period boundary conditions. This is due to the 



absence of any other wave defect, which would eventually destroy 
the source. However static, the source still exhibits interesting dy- 
namics. The wavelength of the waves it generates can vary, as at 
t = 1400 orbits when the wavelength drops abruptly to a much 
shorter value. The cause of this behaviour is unclear, as is its con- 
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nection to the preferred wavelength of the system, which we have 
argued earlier is set by the first linearly stable wavetrain. 

We believe that such a setup might be the most realistic simu- 
lation of the long term behaviour of the overstability. No spurious 
effect from the periodicity can dominate the nonlinear dynamics. 
Note that in this setup it is the shock/source state that is favoured; 
not the earlier disordered wave turbulence. It is then most likely 
that it is the source/shock state that dominates in the real rings. 

4.3.2 Surface density gradients 

We know that the surface density in Saturn's rings is not constant. 
For example, density waves can create enhancements and deficits 
in the local optical depth that might affect the overstability in a 
realistic ring. 

To approximately account for these features, we ran simula- 
tions with an initial linear surface density gradient varying from 
0.04 to 3.24. At the right edge, we sharply drop the optical depth to 
zero. Therefore the right edge is again expected to viscously spread 
(see previous section). We present a space-time plot of our run in 
Figure [9(b)l 

As in earlier runs, the overstability reaches nonlinear ampli- 
tudes within a few hundred orbits. Note that the ring is initially sta- 
ble where x < -7000 m because the optical depth is too low. Sim- 
ilar to the open boundary case presented above, one source region 
eventually wins and is maintained for the rest of the simulation. 
The waves that it generates march leftwards into the lower optical 
depth region and steadily decrease in amplitude and wavelength, an 
interesting effect analogous to that shown in Section |4~2] Note that 
the waves seem to bend upwards on the left side of the plot. This is 
because the wavelength/group velocity depends on optical depth. 

4.4 Synthetic observations 

Although this paper is focused on the numerical and theoretical 
study of the viscous overstability, we present one result from a syn- 
thetic stellar occultation observation. In such an observation, the 
spacecraft looks through the rings onto a background star. The star 
has a finite size and because both Saturn and the spacecraft move, 
the star effectively moves across the ring. By measuring the photo- 
metric optical depth f or the transmission (the fraction of the light 
not blocked by ring particles) one can estimate the local optical 
depth. 

Let us assume an effective stellar radius of r t = 70 m, an inte- 
gration time of the photometer of At = 60 ms, a sampling length of 
Al = 302 m (this could also be expressed in terms of a relative ve- 
locity) and a path of the star which has a projected angle of tp = 5.7° 
relative to the azimuthal direction. These parameters were chosen 
to allow a comparison with observations that can be made with the 
Cassini spacecraft (M. Hedman, private communication). 

The result of applying this convolution is shown in Figure [TO] 
We use a snapshot of the fiducial simulation after t = 4000 orbits 
and plot f (1- transmission) as a function of the time of measure- 
ment and radial distance (relative to the starting point). The convo- 
lution (integration) using the above parameters makes it different 
from the optical depth plotted in Figures [5(c)| and [5(a)| Sharp fea- 
tures are gone, but one can still easily see the radial structure cre- 
ated by the viscous overstability. The amplitudes are of order ~ 0.4. 
The snapshot was chosen so that not all wave-trains have merged 
and a sources is present just right of the centre in the plot. 



Early data from the Cassini spacecraft suggests that structures 
like these are indeed observed. If confirmed, a more detailed and 
systematic comparison is needed to fully uncover the underlying 
ring structure. This work should be seen as a proof of concept. We 
make the code to produce both the simulations and the synthetic 
observations freely available (see discussion section) so as to jump- 
start future investigations. 



5 DISCUSSION 

We have undertaken the first local A'-body simulations of the vis- 
cous overstability in planetary rings over the large spatial and tem- 
poral scales necessary to fully describe its nonlinear saturation. 
With radial domains extending to 50 km and timespans to 10 4 or- 
bits, our particle simulations of overstability are the largest yet at- 
tempted and can begin to compete with one-dimensional hydrody- 
namical calculations. These great scales are achieved by severely 
narrowing the azimuthal box size, a permitted saving because the 
collective dynamics is essentially axisymmetric, and by new nu- 
merical methods implemented in the REBOUND code. 

The simulated dynamics exhibit a wealth of interesting fea- 
tures on different lengths and times ranging over the dimensions 
of our simulations. In particular, these features appear to corre- 
spond to both the periodic microstructure (~200 m) and irregular 
intermediate variations (~ 1-10 km) observed in Saturn's A and fi- 
rings ( |Porco et al.|2005]|Colwell et al.|2007[|Thomson et al.|2007| >. 
They also can be directly compared to those appearing in large- 
scale hydrodynamical simulations. We may then differentiate ro- 
bust dynamical behaviour common to both and which we may then 
expect to prevail in the real rings. 

The most robust dynamical features in our numerical work are 
axisymmetric nonlinear density wavetrains, propagating either in- 
ward or outward. Their wavelength varies between roughly 200 m 
and 400 m, with the longer waves supported by disks with higher 
optical depth f . We believe that the nonlinear wavetrains can be 
unambiguously connected to the observed periodic microstructure, 
though we recognise that our simulated waves are slightly longer 
that those observed via Cassini 's occultation experiments. In addi- 
tion, we find that multiple wavetrains interact in complicated ways. 
Earlier in the evolution, superpositions of multiple wavetrains form 
nonlinear standing modes and beating patterns ('wave turbulence'). 
On intermediate times, however, these resolve into a cleaner struc- 
ture comprising radial patches of waves separated by wave-defects, 
i.e. source and shock regions. Shock regions are where counter- 
propagating wavetrains collide, sources where they are generated 
and move apart from each other. Each wavetrain patch extends over 
many kilometres, whereas the wave defects are shorter — typically 
about 1 km in radial size. It is possible that this longer structure, 
associated with the sources, shocks, and wavetrain patches, corre- 
sponds to irregular intermediate scale variations observed by the 
Cassini cameras. Note that we do not find structure generated on 
scales longer than 10 km, nor do our longest features resemble the 
observations of 100 km scale features in the B-ring. We conclude 
that these are not generated by viscous overstability. The origin of 
100 km structure may lie in ballistic transport instability, electro- 
magnetic instability, or probably both (Durisen 1995; Latt er et aT] 
[20T2l|Goertz & Morfill|1988l >. 

All our simulations with shear periodic boundaries ultimately 
relax to a similar long term state: a single wavelength wavetrain 
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Figure 10. Synthetic stellar occultation observation as it may be seen by Cassini using the fiducial run after t = 4000 orbits. 



propagating in one direction. This state should correspond to the 
shortest stable wavetrain permitted in the box. However, a single 
uniform wavetrain solution is not a realistic outcome in the real 
rings, and is surely an artefact of the periodic boundaries and the 
neglect of any radial structure, i.e. we assume translational sym- 
metry in radius. We conducted additional simulations that broke 
this symmetry. First we examined a spreading ring, whose edges 
were far from the box's radial boundaries. This set-up approximates 
outflow boundary conditions or buffered boundaries because waves 
that impact on the ring edges do not reflect and are lost. We also 
simulated a disk with radial density gradient. In both cases the fi- 
nal single wavetrain solution was averted, with the system settling 
down to a state characterised by a single source region radiating 
density waves. These simulations suggest that the overstability sat- 
urates in the real rings through an arrangement of couterpropagat- 
ing wavetrains interspersed by sources and shocks, embedded in 
the underlying disk structure. However, there may still exist ring 
regions stuck in wave turbulence, exhibiting different photometric 
properties. 

Overall the results from the particle code are in good qual- 
itative agreement with comparable hydrodynamical simulations 
( |Latter &~O gilvie 2010). The dominance of the traveling wave- 
trains, the emergence of sources and shocks, and the long-term 
single-wavetrain states are dynamical behaviours shared by both 
approaches. Nevertheless, interesting discrepancies exist. Impor- 
tantly, particle simulations of colliding wavetrains show greater in- 
terpenetration of the waves than in hydrodynamics where, in gen- 
eral, colliding waves are destroyed immediately. This leads to much 
broader shock regions and a longer and more complicated early 
evolution involving standing waves. Source regions are also wider, 
while the low level chaotic variability in amplitude and phase wit- 
nessed in hydro is much exacerbated in the A'-body dynamics. 

Lastly, we performed a brief set of synthetic occultation obser- 
vations of our numerical output. Direct comparison of such results 
with real occultation data (such as from Cassini UVIS) could in 
principle allow the measurement of a range of particle properties, 
such as size, density, and e. One extension of this work could in- 
volve synthetic observations of reflected light, such as those calcu- 
lated by[S3o"&Kaij3mneii] ^2003j>, |SaIoetd]([2004|, and |Salo &| 
|French|(2010) who successfully explain observations of azimuthal 
brightness asymmetry, the opposition effect, and the tilt effect. Ap- 
plying these techniques to the simulated intermediate scale struc- 
ture would directly connect overstability to the observations on 1- 
10 km from Cassini Imaging Science ( [Porco et al.|2005| l. 

Future work will also assess the influence of self-gravity. First, 



how does it change the axisymmetric dynamics, such as the proper- 
ties of the wavetrains and the larger-scale structure that play upon 
them? Second, how does it change the non-axisymmetric dynamics 
via self-gravity wakes? On larger domains and over longer times, 
can the wakes suppress overstability or vice-versa? What is the na- 
ture of their interference? These last questions also lead to issues 
concerning the longitudinal coherence of the collective axisymmet- 
ric dynamics. Over what azimuthal extent do the density waves 
(and other features) remain coherent? What is the characteristic co- 
herence length and what determines it? Obviously, to answer these 
questions we must move on from the small L v regime employed 
here. 

We make the entire source code publicly available under the 
GPL license. We want to encourage other scientists to both re- 
produce and expand upon our results. The source files are pro- 
vided as a git repository and hosted on github at http : / /github . 
com/hannorein/rebound/tree/overstability See |Rein &| 
|Liu| ( |20i2| > on further informations about the REBOUND code. Fur- 
thermore, the authors are happy to answer any questions that might 
arise with the installation. 
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APPENDIX A: EXTREMELY LARGE RADIAL DOMAIN 

In this appendix we present one simulation with an extremely large 
radial domain, spanning over 55 000 particle radii. With the excep- 
tion of the box size, the simulation parameters are identical to the 
fiducial simulation (see also Tab.[TJ. A stroboscopic space-time di- 
agram of the first 6500 orbits is presented in Fig. |A1| Once the 
overstability starts growing three source regions can be identified 
initially at / ~ 1000 orbits. These are located at x = -20000 m, 
x = 5 000 m and x = 18 000 m. A few thousand orbits later, 
t ~ 4000 orbits, only two survive. At the end of the simulation at 

t > 5000 orbits, only one source can be seen near x 26 000 m. 

This behaviour is qualitatively identical to the fiducial simulation 
presented in Sect. |4.1| 
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Figure Al. Simulation with an extremely large radial domain. 
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